library(ggplot2)
library(patchwork)

# set to desired working directory
# setwd("")

# Start log file for script -----------------------------------------------
TeachingDemos::txtStart("./log files/figureA21.txt")

# Load Data ---------------------------------------------------------
# individual data
load("./data/impdat20.Race.WeekWeight.rdata")
load("./data/impdat20.Age.WeekWeight.rdata")
load("./data/impdat20.Arab.WeekWeight.rdata")
load("./data/impdat20.Asian.WeekWeight.rdata")
load("./data/impdat20.Disab.WeekWeight.rdata")
load("./data/impdat20.GendCar.WeekWeight.rdata")
load("./data/impdat20.GendSci.WeekWeight.rdata")
load("./data/impdat20.NatAm.WeekWeight.rdata")
load("./data/impdat20.Pres.WeekWeight.rdata")
load("./data/impdat20.Religion.WeekWeight.rdata")
load("./data/impdat20.Sexuality.WeekWeight.rdata")
load("./data/impdat20.Weapons.WeekWeight.rdata")
load("./data/impdat20.Weight.WeekWeight.rdata")

# * Age ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Young_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Age.WeekWeight, 
                    subset =  date_mdy != "2020-05-25",
                    weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Age"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"

# Plotting
age_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Arab ------------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.OtherPeople_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Arab.WeekWeight, 
                    subset =  date_mdy != "2020-05-25", 
                    weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Arab"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
arab_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Asian -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.White_American_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Asian.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Asian"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")


preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
asian_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Disability --------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Abled_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Disab.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Disability"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
disability_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())



# * Gender-Career ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Male_Career_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.GendCar.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Gender-Career"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
gendcar_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Gender-Science ----------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Male_Science_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.GendSci.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Gender-Science"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")


preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
gendsci_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Native American ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.White_American_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.NatAm.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Native American"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
natam_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# President ---------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Trump_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Pres.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "President"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")


preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
pres_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Religion ----------------------------------------------------------------

# ^ Christianity-Islam ----------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Christianity_Ci_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Christianity-Islam"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
relig_ci_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# ^ Christianity-Judaism --------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Christianity_CJ_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Christianity-Judaism"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")


preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
relig_cj_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# ^ Judaism-Islam ---------------------------------------------------------

# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Judaism_Ji_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Judaism-Islam"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"

# Plotting
relig_ji_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Sexuality ---------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Straight_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Sexuality.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Sexuality"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")


preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"


# Plotting
sexuality_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Weapons -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Black_Weapons_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Weapons.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Weapons"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"

# Plotting
weapons_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Weight -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_full_dummy <- lm(D_biep.Thin_Good_all ~ 
                      as.factor(week) + ideo_mod + ideo_con + 
                      age_cat4 + edu_cat4 + race5 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Weight.WeekWeight, subset =  date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_full_dummy)

# ** Figure ----------------------------------------------------------
title <- "Weight"
weeks <- names(table(m1_full_dummy$model$`as.factor(week)`))

pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds <- predict(m1_full_dummy, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"

# Plotting
weight_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Race IAT for Placebo Comparison -----------------------------------------
m1_full <- lm(D_biep.White_Good_all ~ 
                as.factor(week) + ideo_mod + ideo_con + 
                age_cat4 + edu_cat4  + sex + 
                census_region + broughtwebsite2,
              data = impdat20.WeekWeight, weights = WeekWeights, 
              subset =  date_mdy != "2020-05-25")

weeks <- names(table(m1_full$model$`as.factor(week)`))

# Data Frames for Predictions
pdat <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")



weeks <- names(table(m1_full$model$`as.factor(week)`))
preds <- predict(m1_full, pdat, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    pred = preds$fit,
                    se = preds$se.fit)
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1)] <- "2020-01-01"

# Plot
race_fig <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2, color = "grey") +
  geom_linerange(aes(ymin = pred - 1.65*se, 
                     ymax = pred + 1.65*se),
                 lwd = .75, color = "grey") +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred), se = FALSE, color = "black") +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred), se = FALSE, color = "black") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "") +
  scale_y_continuous(labels = numform::ff_num(zero = 0, digits = 2)) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 20),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 20),
        strip.background = element_blank())

lhs_fit <- ggplot_build(race_fig)$data[[3]]
lhs_fit$y[which(lhs_fit$PANEL == 1)]

rhs_fit <- ggplot_build(race_fig)$data[[4]]
rhs_fit$y[which(rhs_fit$PANEL == 1)]


race_preds1 <- data.frame(iat = "race",
                          out = "loess",
                          diff = rhs_fit$y[which(rhs_fit$PANEL == 1)][1] - lhs_fit$y[which(lhs_fit$PANEL == 1)][length(lhs_fit$y[which(lhs_fit$PANEL == 1)])])
race_preds2 <- data.frame(iat = "race",
                          out = "predicted",
                          diff = p_dat[which(p_dat$week == "2020-05-25"), "pred"] - p_dat[which(p_dat$week == "2020-05-18" ), "pred"])
race_preds <- rbind(race_preds1, race_preds2)
race_preds <- reshape::melt(race_preds, id = c("iat", "out"), variable_name = "diff")


# Plot -----------------------------------------------------
weeks <- unique(p_dat$week)

obs_diffs <- as.data.frame(array(NA, c(14, 2)))
names(obs_diffs) <- c("iat", "diff")

placeb_diffs <- as.data.frame(array(NA, c(1, 3)))
names(placeb_diffs) <- c("iat", "week", "diff")

x_diffs <- as.data.frame(array(NA, c(length(weeks), 3)))
names(x_diffs) <- c("iat", "week", "diff")

plot_names <- grep("_plot", ls(), value = T)
plot_names <- plot_names[which(plot_names != "facet_plot")]

# loop across plot objects to pull predicted D-scores for observed differences and pseudo differences
for(i in 1:length(plot_names)){
  # metadata storing
  x_diffs$iat[1] <- gsub("_plot", "", plot_names[i])
  x_diffs$week[1] <- weeks[1]
  
  # pull observed George Floyd Difference
  obs_diffs$iat[i] <- gsub("_plot", "", plot_names[i])
  
  lhs_fit <- ggplot_build(get(plot_names[i]))$data[[1]]
  rhs_fit <- ggplot_build(get(plot_names[i]))$data[[2]]
  
  obs_diffs$diff[i] <- rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == "2020-05-25")[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == "2020-05-18")[1]]

  # other differences but for the treatment period
  for(j in 2:length(weeks)){
    x_diffs$iat[j] <- gsub("_plot", "", plot_names[i])
    x_diffs$week[j] <- weeks[j]
    
    if(weeks[j] <= "2020-05-18"){
      x_diffs$diff[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == weeks[j-1])[1]]
    }
    if(weeks[j] == "2020-05-25" | weeks[j] == "2020-06-01"){
      # empty placeholder
    }
    if(weeks[j] > "2020-06-01"){
      x_diffs$diff[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == weeks[j-1])[1]]
    }
  }
  
  # store results
  placeb_diffs <- rbind(placeb_diffs, x_diffs)
}

# Reformatting and cleaning up
placeb_diffs <- placeb_diffs[-1,]
placeb_diffs <- reshape::melt(placeb_diffs, id = c("iat", "week"), variable_name = "diff")
obs_diffs <- reshape::melt(obs_diffs, id = c("iat"), variable_name = "diff")

# Getting p-values 

# Data frame to store p-values
p_values <- obs_diffs

for(i in 1:nrow(p_values)){
  # store test statistic
  ob_val <- obs_diffs$value[i]
  # extract appropriate simulations for reference distribution
  ref_dist <- subset(placeb_diffs, iat == p_values$iat[i],
                     select = value)
  ref_dist <- ref_dist$value[which(!is.na(ref_dist$value))]
  # store number of simulations
  n_sim <- length(ref_dist)
  # two-tailed p-value
  p_values$value[i] <- sum(abs(ref_dist) > abs(ob_val))/n_sim
}
p_values$x <- -.05
p_values$y <- 10
p_values$p <- paste0("p=", round(p_values$value, 2))

unique(placeb_diffs$iat)

labs <- c("Age", "Arab", "Asian", "Disability",
          "Gender-Career", "Gender-Science",
          "Native\nAmerican", "President", 
          "Religion\nChristian-Islam", "Religion\nChristian-Judaism",
          "Religion\nJudaism-Islam", "Sexuality", "Weapons", "Weight")

placeb_diffs$iat <- factor(placeb_diffs$iat, labels = labs)
p_values$iat <- factor(p_values$iat, labels = labs)
obs_diffs$iat <- factor(obs_diffs$iat, labels = labs)


# Plotting
ggplot(placeb_diffs, aes(x = value)) +
  geom_histogram(fill = "grey") +
  geom_vline(data = obs_diffs, aes(xintercept = value)) +
  geom_text(data = p_values, aes(label = p, x = x, y = y)) +
  facet_wrap(iat~.) +
  scale_x_continuous(labels = numform::ff_num(zero = 0, digits = 2), limits = c(-.1, .1)) +
  theme_bw() +
  labs(x = "Predicted IAT Score Difference", y = "Count",
       # title = "Comparing Predicted Pre-Post Race IAT Change to Weekly Variation in D-Scores by IAT Type",
       # subtitle = "Observed Difference is Week of May 18 to Week of May 25"
  ) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 16),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 14),
        strip.background = element_blank())
ggsave("./figures/FigureA21.pdf", height = 8, width = 8)

# End log file for script -------------------------------------------------
TeachingDemos::txtStop()